Feedback control of heart rate during treadmill exercise based on a two-phase response model

This work investigated automatic control of heart rate during treadmill exercise. The aim was to theoretically derive a generic feedback design strategy that achieves a constant input sensitivity function for linear, time-invariant plant models, and to empirically test whether a compensator C2 based on a second-order model is more dynamic and has better tracking accuracy than a compensator C1 based on a first-order model. Twenty-three healthy participants were tested using first and second order compensators, C1 and C2, respectively, during 35-minute bouts of constant heart rate treadmill running. It was found that compensator C2 was significantly more accurate, i.e. it had 7% lower mean root-mean-square tracking error (1.98 vs. 2.13 beats per minute, p = 0.026), and significantly more dynamic, i.e. it had 17% higher mean average control signal power (23.4 × 10−4 m2/s2 vs. 20.0 × 10−4 m2/s2, p = 0.011), than C1. This improvement likely stems from the substantially and significantly better fidelity of second-order models, compared to first order models, in line with classical descriptions of the different phases of the cardiac response to exercise. These outcomes, achieved using a treadmill, are consistent with previous observations for the cycle ergometer exercise modality. In summary, whenever heart rate tracking accuracy is of primary importance and a more dynamic control signal is acceptable, the use of a compensator based on a second-order nominal model is recommended.


Introduction
During exercise training, it can be beneficial to vary exercise intensity between two or more levels for varying durations: this so-called "interval training" has advantages compared to continuous training [1][2][3].Because heart rate (HR) is a variable that is commonly used to set exercise intensity [4], this has inspired the development of accurate and robust HR control systems.In the present work, we employed a treadmill (TM) as the platform for HR control.Compared to other exercise modalities, e.g.cycle ergometry (CE), treadmill exercise has a relatively high energy expenditure at a given HR [5], especially at low to moderate intensities, cf.[6].
Recent studies of feedback control of heart rate are mostly model-based and can be divided into two classes: several investigations [7][8][9][10] employed the nonlinear model proposed by Cheng et al. [11], while others did the feedback design using linear time-invariant models [12,13].Because the stability of nonlinear controllers is conditional [8,10], the corresponding HR control studies mainly focused on theoretical stability and robustness analysis, while empirical verification used either simulation [7,8] or experiments with very small numbers of participants [9,10].
In contrast, recent studies with large and statistically meaningful test group sizes showed that a first-order, linear, time-invariant (LTI) model can be used to design LTI compensators that achieve impressive HR control performance: a HR control study on healthy participants (n = 30) achieved a mean root-mean-square heart rate tracking error (RMSE) of 2.96 beats per minute (bpm) [14]; another study on healthy participants (n = 25) that compared linear HR control performance on TM and CE modalities achieved a mean RMSE of 2.85 bpm on the treadmill and 3.10 bpm on the cycle ergometer [15].Moreover, a systematic comparative study with sample size n = 16 showed that HR control performance with linear and nonlinear controllers was not significantly different; in fact, the nonlinear controller had worse performance at low speeds [16].
With a view to further improving the performance of HR control based on linear models, a second-order model structure was investigated in our previous treadmill-based model identification study [17].This structure was motivated by physiological knowledge regarding human HR response to exercise that identifies three phases: Phase I is a small but immediate response, Phase II is a large and slower component that comprises the main part of the overall response, and Phase III is an ultra-slow drift that can occur above the anaerobic threshold [18].Since, in feedback control, very slow disturbances can be effectively eliminated by integral action, Phases I and II are significant for feedback design: this motivated our investigation of twophase, i.e. second order, models.The comparative analysis with n = 22 showed that the second-order model structure achieved significantly better fit due to the explicit delineation of the fast and slow dynamic components of HR response.Furthermore, the improvement is not limited to treadmills: a similar model identification study using a cycle ergometer (n = 26) gave consistent results [19].
The improved fidelity of second-order models, due to the inclusion of fast and slow dynamics, leads to the hypothesis that compensators based on these models could be more dynamic and thus might deliver more accurate heart rate tracking.The concept of a controller being more or less dynamic is quantified in this work using the average power of changes in the control signal (the treadmill speed command; see Eq (32) in the sequel).This hypothesis was experimentally tested in a pilot study with 10 participants [20].It was found that compensators derived using second-order models, denoted C 2 , were significantly more dynamic than compensators based on first-order models, C 1 .However, there were no significant differences in heart rate tracking accuracy between C 1 and C 2 .This outcome was likely due to important differences between the sensitivity functions for the two compensators, differences in the reference prefilter design, and the sample size being too small to allow detection of differences, even when they exist.These limitations were addressed in a subsequent study using a cycle ergometer [19]: with n = 26 participants, controller type C 2 was found to be significantly more dynamic and more accurate than type C 1 .
This analysis motivates further study of control performance of C 2 vs. C 1 using the treadmill modality.Here, we rectify the limitations of the treadmill pilot study [20], in the following way: 1. the feedback design was modified to make the input-sensitivity function constant over all frequencies, which in consequence made the other sensitivity functions flatter (less peaking) and the nominal functions for C 1 and C 2 more similar; this required a new theoretical derivation, presented here, for the second-order case of a feedback compensator C 2 that gives a constant input sensitivity function; 2. the confounding effect of the reference prefilter was eliminated by considering only a constant HR reference signal; 3. statistical power was improved by increasing the sample size to be similar to previous studies that were well powered, viz.we used n = 23.
The aim of this work was twofold: (i) theoretical contribution-to develop a novel feedback design strategy that achieves a constant input sensitivity function for LTI plant models in general, and for the second-order case in particular; (ii) experimental contribution-to empirically test whether a compensator based on a second-order model is more dynamic and has better HR tracking accuracy than a compensator based on a first-order model.

Controller design
The control structure used for this study has a standard form consisting of two parts (Fig 1 ): a feedback compensator C(s) that adjusts the control signal u, i.e. the commanded speed of the treadmill, in dependence upon the error e between filtered HR reference r 0 and the HR measurement z; and a reference prefilter C pf to manipulate the overall tracking response of the closed-loop system, i.e. the response from reference heart rate r = HR* to actual heart rate y = HR.A disturbance term d models heart rate variability (HRV) and other sources of uncertainty, and n represents measurement noise.
Nominal plant model.As detailed in our previous model identification study [17], the general form for the nominal plant P o is taken to be the strictly proper transfer function In the sequel we solve the general compensator problem for the generic plant Eq (1) and then specialise the solution to two instances of the plant: a first-order form where B(s) = k 1 /τ 1 and A(s) = s + 1/τ 1 , and a second-order form where where the G, H and H 0 have degrees n g , n h and n h 0 , respectively, and H is taken to be monic.The integral term 1/s is introduced to compensate the very slow Phase III dynamic of HR response.
In the following, we consider the transfer function C to be merely proper (rather than strictly proper) with n g = n h .This is to prevent the input sensitivity function U o (Eq (7), below) from rolling off to zero with increasing frequency, in line with our design goal to maintain a constant |U o (jω)| for all frequencies: since Using the rational forms for the plant P o , Eq (1), and compensator C, Eq (4), the input-sensitivity function can be written where F = AH + BG is the closed-loop characteristic polynomial.
Compensator synthesis for constant input sensitivity: General solution.Since the feedback design goal is to achieve an input-sensitivity magnitude that is constant for all frequencies, it follows that the closed-loop poles (the roots of F) must be completely cancelled by the numerator polynomial AG in Eq (8).This can be achieved by firstly cancelling the plant poles by including A in the compensator numerator, i.e. by setting G = AG 0 , giving Second, the remaining closed-loop poles (the roots of H + BG 0 ) are also placed at the openloop poles by setting H + BG 0 = A to give U o = G 0 .Considering now that we have set F = A 2 with degree n φ = 2n a , and that the numerator AG = A 2 G 0 of Eq (8) must have this same degree, it follows that 2n a þ n 0 g ¼ 2n a , thus n 0 g ¼ 0 and G 0 is therefore a constant, denoted g 0 0 .Consequently, and which, as desired, is a constant value.The compensator parameters are obtained by solution of Eq (10) for g 0 0 and H 0 ; since the compensator is proper (n g = n h ), and since G ¼ g 0 0 A, n g = n h = n a and The controller solution derived above is valid for the generic nominal plant model P o = B/A, Eq (1), without any restriction on the plant order.In summary, with the constraint G ¼ g 0 0 A, the compensator Eq (4) that achieves a constant U o ¼ g 0 0 , Eq (11), is given by where constant g 0 0 and polynomial H 0 are the unique solution of the polynomial Eq (10).The cancellation strategy also leads to simplifications in the expressions for S o and T o .From Eqs ( 5) and ( 6), and using Compensator synthesis for constant input sensitivity: First and second order solutions.In the following, the general solution is specialised to first and second-order plants described by Eqs (2) and (3), respectively.
First-order solution.In the first-order case with n a = 1, Eq (12) gives n 0 h ¼ 0 which in turn leads to the trivial solution H 0 = 1 and H = s (H, and therefore also H 0 , were assumed at the outset to be monic).Eq (10) then reads which, by equating coefficients, has the solution For first-order plants, the compensator is therefore (see Eq (13)) and the corresponding sensitivity functions are, using Eqs ( 11), ( 14) and ( 15), Second-order solution.In the second-order case with n a = 2, Eq (12) gives A solution is then sought for Eq (10), which is obtained by straightforward algebraic manipulation as For second-order plants, the compensator is therefore (see Eq (13)) From Eqs (11), ( 14) and ( 15), the corresponding sensitivity functions are The Eqs ( 18) and ( 24) for C 1 and C 2 show that the compensators comprise very simple expressions that depend only on the nominal plant parameters: no additional tuning parameters are otherwise involved.
Compensator calculation.For calculation of the compensator coefficients, nominal first and second-order plant parameters obtained from our previous identification study [17], were used: these models are averages obtained from 22 individual first and second order models from 11 participants (for each participant there were two first-order models and two secondorder models), with values k 1 = 28.57bpm/(m/s), τ 1 = 70.56s and k 2 = 24.70 bpm/(m/s), τ 21 = 18.60 s, τ 22 = 37.95 s.Substituting in Eqs ( 18) and ( 24), the compensators are and The magnitude plots of the resulting closed-loop sensitivity functions for the pairs P 1 , C 1 and P 2 , C 2 are illustrated in Fig 2 .It can be seen that the input sensitivity magnitudes are constant, according to Eqs ( 19) and (25), and that the complementary sensitivity function for the second-order case rolls off twice as fast as for the first-order case (40 dB/decade vs. 20 dB/ decade, Eq (21) vs. Eq (27)).
Finally, the reference prefilter C pf was designed to make the overall closed-loop transfer function (r 7 !y, see Fig 1) equal to a second-order transfer function denoted where T o is the nominal complementary sensitivity function in each case.T o had critical damping and 150 s rise time (detailed in [14]).The transfer function of C pf was then o ðsÞ � T cl ðsÞ : r 7 !r 0: ð30Þ

Experimental design
Twenty-three healthy participants were recruited for this study with ages between 23 years and 57 years, mean body mass 71.4 kg and mean height 1.77 m.All participants were regular exercisers (at least 3 times a week, 30 minutes each time), non smokers and free from cardiovascular disease or musculoskeletal complaints.The study was approved by the Ethics Committee of the Swiss Canton of Bern (Ref.2019-02184) and participants provided written, informed consent prior to participation.Each participant completed two feedback control tests on the treadmill, one with C 1 and one with C 2 .The order of presentation of each compensator (C 1 first or C 2 first) was changed in sequence for each participant to counterbalance possible order-of-presentation effects.Each test was conducted on a separate day and consisted of a 10-minute warm up, 10-minutes rest and a 35-minute formal measurement (Fig 3).During the warm up and formal measurement phases, the treadmill speed was controlled by the HR compensator (C 1 or C 2 ) to regulate the HR of the participant at a constant, mid-level target value, denoted HR m , that was set to an intensity level at the border between moderate and vigorous using the prediction equation HR m = 0.765 × (220 − age) [23].The 30-minute time interval from 295 s to 2095 s was used for outcome evaluation ("evaluation period" in Fig 3), except for participants P04 and P09 where the period 295 s to 1795 s was used (for both of these participants, there were measurement artefacts in the final five minutes of the tests).

Equipment
We used a PC-controlled treadmill (model pulsar, h/p/cosmos Sports & Medical GmbH, Germany; Fig 4).The heart rate compensators were implemented in Simulink Desktop Real-Time Closed-loop frequency responses.Input-sensitivity functions U 1,2 , sensitivity functions S 1,2 and complementary sensitivity functions T 1,2 where the 1 and 2 subscripts denote the first and second-order plant/ compensator combinations P 1 , C 1 and P 2 , C 2 , respectively.The four vertical dashed lines at 0.00 � 3 Hz, 0.04 Hz, 0.15 Hz and 0.4 Hz delineate the four frequency bands classically used for heart rate variability analysis [22]: ultra-low frequency (ULF), very-low frequency (VLF), low frequency (LF) and high frequency (HF).https://doi.org/10.1371/journal.pone.0292310.g002(The MathWorks, Inc., USA) and speed commands were sent to the treadmill over a serial communication protocol.Heart rate was measured by a chest belt sensor (H10, Polar Electro Oy, Finland) and sent to Simulink through Bluetooth.Heart rate was sampled at 1 Hz while the compensators ran at a rate of 0.2 Hz (sample interval 5 s), hence the heart rate was downsampled by averaging every five consecutive values.

Outcome measures and statistical analysis
The accuracy of heart rate tracking was quantified using the root-mean-square error (RMSE) between measured HR and the nominal, simulated HR response.The intensity of the control signal, i.e. the treadmill speed command u, was evaluated by the average power of changes in this variable (P ru ): RMSE ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 N X N i¼1 ðHR sim ðiÞ À HRðiÞÞ HR sim represents the simulated closed-loop HR response, i.e.HR sim = T cl � HR*, but since the target HR was constant during these experiments, HR sim = HR m during the evaluation period: that is to say, the reference prefilter only played a role during the initial transient prior to the evaluation period and had no effect on the outcome calculation.i indexes the discrete sample instants and N is the total number of samples during the evaluation period.As the evaluation period was from 295 s to 2095 s (Fig 3) and the sample rate of the compensator was 0.2 Hz, N = 361 (except for P04 and P09, where N = 301 as the evaluation period was from 295 s to 1795 s).
The hypothesis of this study was that, compared with C 1 , C 2 would have better HR tracking accuracy, signified through a lower RMSE, and a more dynamic control signal, manifest by a higher P ru ; thus, for both outcomes, one-sided hypothesis testing was done.Prior to hypothesis testing, a Kolmogorov-Smirnov test with Lilliefors correction was used to assess normality of sample differences.As all the differences were found not to significantly deviate from normality, testing was performed using paired t-tests.The significance level was set as α = 0.05.Statistical calculations were implemented in the Matlab Statistics and Machine Learning Toolbox (The Mathworks, Inc., USA).

Discussion
The aim of this work was to theoretically derive a generic feedback design strategy that achieves a constant input sensitivity function for LTI plant models, and to empirically test whether a compensator based on a second-order model is more dynamic (higher P ru ) and has better HR tracking accuracy (lower RMSE) than a compensator based on a first-order model.
The theory for constant input sensitivity feedback design was first derived for general LTI plant models and was then specialised to the first and second-order cases.For the specific compensators employed here, it was found that the input sensitivity functions U o were indeed constant for all frequencies, and that the sensitivity and complementary sensitivity functions S o  18) and ( 24), have a very simple functional form and their coefficients depend only on the respective nominal model parameters (gains and time constants).
Overall, it was found that compensator C 2 was significantly more accurate, i.e. it had 7% lower mean RMSE, and significantly more dynamic, i.e. it had 17% higher mean P ru , than C 1 .This improvement likely stems from the substantially and significantly better fidelity of second-order models compared to first order models [17], in line with classical descriptions of the different phases of the cardiac response to exercise [18].These outcomes, achieved using a treadmill, are consistent with previous observations for the cycle ergometer exercise modality [19].
In summary, whenever heart rate tracking accuracy is of primary importance and a more dynamic control signal is acceptable, the use of a compensator based on a second-order nominal model is recommended.

Fig 1 .
Fig 1.Control structure of this study.Nominal plant model P o (s), feedback compensator C(s) and reference prefilter C pf (s).The controlled variable y is the heart rate (HR) and r is the reference/target heart rate (HR*).The control signal u is the commanded speed of the treadmill.The terms d and n respectively represent disturbances and measurement noise.https://doi.org/10.1371/journal.pone.0292310.g001 and since P o is by definition strictly proper (whence lim ω!1 |P o (jω)| = 0), U o � C at high frequency.With a proper C, U o must also be proper and |U o (jω)| therefore remains finite.For the feedback system of Fig 1, the principal sensitivity functions-sensitivity function S o , complementary sensitivity function T o and input-sensitivity function U o -are defined classically [21], as

Fig 5 .
Fig 5. Sample measurements with C 1 and C 2 that had the lowest, middle and highest RMSE.The upper plot of each figure shows the reference HR (HR*, black dashed line), the measured HR (HR, red line) and the simulated HR response (HR sim , black line).The lower plots show the treadmill speed command (black line).The green horizontal line marks the nominal evaluation period for outcome calculation.Pxy denotes the participant number.A: Lowest RMSE for C 1 , P07.B: Lowest RMSE for C 2 , P12.C: Middle RMSE for C 1 , P22.D: Middle RMSE for C 2 , P14.E: Highest RMSE for C 1 , P10.F: Highest RMSE for C 2 , P06.

Graphical depiction of statistical comparison between outcomes
. Blue (C 1 ) and red (C 2 ) dots are the individual samples; green lines connect samples from individual participants; red bars mark the sample means (given numerically in Table1).D shows the differences between paired samples (C 2 -C 1 ).MD shows mean differences (red horizontal bars) and the corresponding 95% confidence intervals (CIs, blue lines and arrows).When the difference is significant, the value 0 lies outwith the CI; the * notation indicates p < 0.05.A: RMSE for C 1 and C 2 , p = 0.026.B: P ru for C 1 and C 2 , p = 0.011.https://doi.org/10.1371/journal.pone.0292310.g006andTodisplayed little or no peaking (Fig2).Furthermore, the respective closed-loop frequency responses for the C 1 , P 1 and C 2 , P 2 combinations were nominally very similar.It is notable that the compensator transfer functions for both C 1 and C 2 , Eqs (